Jet reconstruction in hadronic collisions by Gaussian filtering 
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A new algorithm for jet finding in hadronic collisions is presented. The algorithm, based on a 
Gaussian filter in (r],4>), is specifically intended for use in heavy ion collisions and/or for detec- 
tors with limited acceptance. The performance of the algorithm is compared to two conventional 
algorithms, a seedless cone algorithm and a k± algorithm, for Pythia simulated di-jet events in 
•Js = 200 GeV p + p collisions with 4GeV/c < a/Q^ < 16GeV/c. The Gaussian filter is found to 
perform as well as, and in some instances better than, the conventional algorithms. 
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I. INTRODUCTION 



The observation of jet quenching in heavy ion colli- 
sions at RHIC (l| and the problems inherent in the use 
of single or di-hadron observables to quantify the en- 
ergy loss of hard-scattered partons [2| together provide 
a strong motivation for complete jet measurements in 
heavy ion collisions. However, at RHIC energies, the 
underlying heavy ion event makes application of conven- 
tional jet algorithms difficult. Commonly used jet algo- 
rithms employ either fixed-size cones [3], augmented by 
split /merging procedures, or iterative clustering 
techniques [J, [5( to determine the solid angle coverage 
of the jets. The kinematic parameters of the jets are 
obtained by summing over the particles or calorime- 
ter elements within the jet with a constant (i.e. angle- 
independent) weight. Such a flat weighting may not be 
optimal in the presence of a fluctuating background be- 
cause a typical jet has most of its energy concentrated 
near the center of the jet and only a small fraction of the 
energy in the periphery of the jet. In contrast, an angular 
weighting that enhances the center of the jet compared 




FIG. 1: A demonstration of the application of the Gaussian 
filter to a Pythia event. Final state particle pr are plotted on 
the bottom Lego plot. The result of the filter is shown with 
the contour plot on the top surface. Red connecting lines 
indicate reconstructed jet positions. 



to the periphery would provide an improved signal to 
background in the measurement of the jet energy. 

Backgrounds from the underlying event in hadron- 
hadron collisions or heavy ion collisions or from pile-up 
at high luminosity colliders can also distort the jet find- 
ing algorithms themselves. High energy particles from 
the background can prevent or modify the convergence of 
the mean-shift iteration in cone algorithms or can mod- 
ify the order of combination of fragments/proto-jets in 
k± algorithms. An algorithm that finds maxima in the 
the angular distribution of fragment (transverse) ener- 
gies using a weighting function or filter that is strongly 
peaked has the potential to find the jet position with less 
sensitivity to the presence of background. Such an al- 
gorithm also has the advantage of reducing the impact 
of the limited acceptance of certain heavy ion detectors 
used at RHIC and the LHC. 

We describe in this paper a new algorithm for find- 
ing jets and extracting jet kinematic quantities that uses 
a linear, Gaussian filter applied in the space of pseudo- 
rapidity and azimuthal angle, (r],4>). Jets are found as 
local maxima of the filter output. Because the filter is 
strongly peaked, the algorithm is expected to improve 
the reconstruction of jet positions and jet energies in the 
presence of background and to reduce the impact of re- 
stricted aperture on jet measurements. This paper takes 
the first step in exploring the basic characteristics of the 
Gaussian filter jet finder by studying its performance on 
Pythia simulated p + p events and comparing the perfor- 
mance to two conventional algorithms, the SISCone Q 
seedless cone algorithm and an implementation [t| of the 
k± algorithm. 

We note that an angular weighting of final state parti- 
cles has been used previously to define energy flow vari- 
ables [H, and the algorithm described here is inspired 
by that work. Though we focus here on jet finding and 
energy estimation, as with energy flow variables, the fil- 
tered event shape contains more information than just 
the locations and energies of reconstructed jets, A pre- 
vious use of low-pass filtering by Donati et al. [9$ should 
also be recognized. However, that work only used the fil- 
ter to bootstrap other clustering algorithms, while in this 
work, filtering is used for complete jet reconstruction. 
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FIG. 2: Filter reconstructed Pythia jet pr distribution func- 
tion for PT,fc_ L ,PT,siscone € {4,8, 16}GeV/c ± 250MeV/c, 
respectively. Jets reconstructed by different algorithms are 
matched within AR = a/At? 2 + A(j> 2 < 0.1. 

II. ALGORITHM 

Instead of the traditional picture of discrete final state 
particles used by traditional jet reconstruction algo- 
rithms, we would like to consider a set of particles (pt,i) 
of generating the event pt density 

PT{ r n,<t>) = ^PT^{ri--q l )8{(j)-(j)i). (1) 

i 

Defining the jet reconstruction over a pr density has the 
ability to compensate for a non-localized background, 
which is useful to apply the jet reconstruction algorithm 
to heavy ion collisions and background in high luminosity 
hadronic colliders. 

The jet reconstruction procedure can be expressed as 
finding the discrete set J of all jets with the transverse 
momentum p?-, pseudo-rapidity r\ and azimuth <fi, 

J = { (pr, r), 4>) | pr ■— Pt(v> 4>) a local maximum}, (2) 

with the filtered pr density pT(r],4>) being the linear- 
circular convolution 



Pt(v,4>) = 




drj'd<j)'pT(v', 4>')KV ~ v\ <t> ~ <t>') 



on the cylindrical (rj, (f) topology. 

When implemented solely with an iterative maximum 
finding, ([2]) would have implicitly the issue of proper ini- 
tialization, not unlike the seeding dilemma of the cone 
algorithm. However, discrete digital filtering provide an 
efficient mean to calculate Pt{t], 4>) f° r a large number of 
sample points. Thus one could find every possible max- 
imum, i.e. in the analogy of the cone algorithm, achiev- 
ing "seedlessnes" by means of sufficiently sample the en- 
tire (rj, 4>) range with a number of seeds N = N n N^ ^> 



FIG. 3: Pythia jet multiplicity of jets with pr > 2GeV/c, 
with PT,tri g = 8GeV/c for the filter, k±, and SISCone algo- 
rithms 

27rA77/i?Q, with Ro w 0.2 being the characteristic separa- 
tion of close jet pairs, and A77 = 77 max — ?7mm the pseudo- 
rapidity range (typically that of the collider experiment 
calorimetry) . 

Therefore we implement ^ by a multistage algorithm, 
consisting of the following steps: 

1. Accumulate a rectangularly binned px density of 
the event. This can be thought as a pt his- 
togram, or in term of the density distribution (|TJ) 
and modulo a constant normalization, to approx- 
imate {f]i,(f>i) by a discretized {f]i,4>i) f° r each i, 
with f) = [(rj- rj min )N v / A77J Ar//N v + (and [-J 
denoting the floor function) being the discretized 
pseudo-rapidity, and analogously for <p. 

2. Apply the discrete realization of the filter on the 
binned density to obtain the initial pr {fi,<j)) ap- 
proximation for the discrete (fj, (f>) bins or pixels. 
This filter could be cither implemented in the (77, <fi) 
position space using e.g. a finite impulse response 
(FIR) or infinite impulse response (IIR) version of 
the filter, or realized in the Fourier space. 

3. The local maxima are localized by comparing the 
filtered px density of each discrete pixel against 
that of the surrounding pixels. 

4. The set of locally maximal pixel centers (77, cf>) is 
used to initialize a suitable local optimization algo- 
rithm, operating on the unbinned pr density fl}, 
to obtain the true (77, <j>) positions. 

A jet definition like @ is inherently collinearly and in- 
frared safe, as pt{Vj 4 1 ) is a non-iteratively defined quan- 
tity not involving any thresholds or cutoffs, and therefore 
like event shape variables such as thrust, is neither sen- 
sitive to infinitely soft radiation, nor collinear splitting. 
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FIG. 4: Triggered Pythia jet spectrum, left with pr.trig = 8GeV/c, comparing different algorithm and size a selection, and 
to the middle comparing filter a = 0.5 and a = 0.7, respectively, against SISCone R = 2 _1,/2 ~ 0.707 and for different 
PT.trig € {4, 8, f 6} GeV/c 



When implemented properly, the maximum finding only 
depends on the collinearly and infrared safe pt(t],4>). 

For the filtering kernel, we would like to propose a bi- 
variate Gaussian distribution function with the normal- 
ization h(0, 0) = 1, 

h( V ,<l>)=e X p[-(rj 2 +4>l I )/(2a 2 )], (3) 

with ar := 2tt\{4> + 7r)/(27r)] — tt (and [•] denoting the 
ceiling function) being the angularly reduced azimuth. 
The parameter a determines the radial scale of the jet re- 
construction, analogous to the fixed-cone radius R. Fig- 
ure Q] illustrates the Gaussian filter applied to a Pythia 
event, and a comparison of the filter reconstructed jet 
Pt distribution functions for jets that the k± algorithm 
and SISCone reconstruct to pt — 4, 8, and 16GeV/c is 
shown in Figure O 

The mean shift iteration (such as the cone algorithm) 
is closely related to a bound maximization with respect 
to the convolution with a "shadow function" [13, . 
This has the consequence that for the pure jet direc- 
tion finding, the filtering based algorithm with subse- 
quent maximization presented here has an exact corre- 
spondence in the picture of a weighted mean shift itera- 
tion, and therefore can be regarded as a generalized form 
of the cone algorithm. The correspondence also leads 
naturally to the definition of jets as the local maxima in 
©. For the cone algorithm, the shadow function kernel 
is h{r],4>) — max[0, 1 — [rj 2 + <f) 2 )/R 2 ]. Comparing this 
with the functional form of © suggest for the limit of 
narrowly focused jet, the angular behavior of the Gaus- 
sian filter and the cone algorithm (without split/merge) 
would correspond with the choice V2a — R. 

There are several reasons to prefer a Gaussian density 
distribution instead of e.g. the parabolic shadow function 
implied by the cone algorithm. A Gaussian density does 
not possess the undesirable feature of maximum creation 



(or in the language of a cone iteration, stable midpoint 
axes). And when a slow varying background is present, 
the foreground to background ratio is accounted for both 
in the reconstructed jet energy and the directional jet 
finding, since the Gaussian weighting is present in both 
Pt(ji,4>) and its gradient. As a rapidly decreasing test 
function, the Gaussian density possesses the mathemati- 
cal property necessary to regulate the singular (77, <f>) di- 
rectional distribution of an outgoing high-py fragment in 
©■ 

Gaussian filtering has been used extensively in fields 
such as computer vision, and several fast approximations 
are known beside the direct and Fourier space convolu- 
tion. We use a partial fraction expanded form of 
the IIR approximation originally described by Young & 
van Vliet [H, [3 [l5[ , which compared to a fast Fourier 
transform (FFT) realization requires no padding along 
the pseudo-rapidity axis, but at least two filtering passes 
to simulate a circular response in azimuth. It should be 
noted that multi-pass filtering does not exactly reproduce 
in the the angularly truncated Gaussian distribution as 
in ([3]) , but rather the sum of its circular frequency com- 
ponents, and the filtering kernel therefore in fact takes 
the form 

hfafl = (2 7 r)- 1 / 2 | ( 7| e -" 2 /(2- 2 ) t9 3(0/2,e- CT2 / 2 ), 

with $3 being the Jacobi theta function [l6(. The dif- 
ference is however negligible for all meaningful choices 
of a <C 7r, especially at the initial, approximate stage of 
filtering. 

As input to the filter we consider final state particles 
within an acceptance of \r\\ < 8. The discrete filtering 
operates on a N v = 640, iV^ = 256 grid. For the con- 
tinuous maximization, we use the Newton optimization 
algorithm with analytically calculated gradient and Hes- 
sian, and for robustness, modified to handle of indefinite 
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FIG. 5: Pythia di-jet event jet opening angle with pT.trig 
8GeV/c for the filter and k± algorithms 



Hessian by spectral decomposition. 



III. SIMULATION AND EVENT SELECTION 

In order to assess the performance of jet reconstruc- 
tion by filtering, we compared the reconstructed jets 
against the k± and cone algorithms at the ideal detec- 
tor level, by using event generator truth particles from 
Pythia 6.4.16 

For all comparisons we present here, we consider the 
p+p collision system at = 200 GeV, which as of Run- 
8 accounts for approximately 170 pb -1 of the hadronic 
RHIC operating mode. The k± algorithm used is the Kt- 
Jet [7j implementation in the inclusive clustering mode, 
with the inclusive stopping parameter D = 1 and the lon- 
gitudinally invariant QCD distance scheme We fur- 
ther compare against SISCone [6j, a variant of the cone 
algorithm, with the overlap threshold parameter / = 0.5 
and the cone radius R = 2~ 1 / 2 0.707, close to the 
typically used R = 0.7 and in asymptotic angular corre- 
spondence to the Gaussian filter with a = 0.5 for narrow 
jets. Each of the comparison presented here are based on 
10 6 (triggered) Pythia events. 

Two form of triggering were used. For the event multi- 
plicity, di-jet and three jet balance (Figures[3j[5l[6j[7]), we 
are comparing two algorithms at a certain scale of hard 
scattering. Here we triggered on the Pythia hard scatter- 
ing Pt.2~>2 with a centered, 500 MeV wide window, i.e. 
withp T , trig -0.25 GeV/c < p T ,2~>2 < p T ,tri g +0.25 GeV/c, 
and obtained the differential cross sections per unit trig- 
ger window at the triggered PT,trig- 

The same applies to 
the triggered pt spectrum that illustrates the sensitivity 
to soft QCD background (Figure [4J . For the compari- 
son of filter reconstructed px against the reconstructed 
Pt of the k± and cone algorithms (Figure ^ , we used a 
large trigger window with pt. 2^2 > PT.tiig ~ 2.25 GeV/c 



FIG. 6: Pythia di-jet event pr sum with PT,tri g = 8 GeV/c for 
the filter and fcj_ algorithms, with the k± convolved by the 
filter to k± energy scale ratio. A Gaussian and exponential 
fit for the low and high inbalance, respectively, are shown. 



to avoid biasing the jet, and the results are shown as 
normalized probability density functions. 

To test the di-jet balance of the filtering algorithm, 
we arrange the reconstructed jets in descending p?, i.e. 
Pt,i > Pt,2 > ■ ■ • > Pt,n- Events with a prominent di- 
jet structure can be selected by requiring pr,3 < 9Pt,2, 
with g < 1 as the relative pr gap, which we chose to 
be j. This approach avoids sensitivity to the individual 
energy scales of the algorithm. In the pairwise compari- 
sion, the selection is applied on both jet algorithms in the 
logical "and" sense to ensure symmetric event selection 
(i.e. both algorithms agree that a clear dijet structure is 
present) . 

An analogous selection for the three jet events can be 
made by requiring pta < 9Pt,3 (with the same g = \ 
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FIG. 7: Pythia three jet event jet opening angle with pr.trig 
8 GeV /c for the filter and k± algorithms 
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choice as with the di-jet event selection). If the pair- 
wise angle among the three leading jets are </>i2, 4>23> </>3i > 
the resulting three jet opening angle is defined as A<fi := 
min(0i2, 4>23, 031 )■ Since three algorithms are compared 
in Figure [71 we used the logical "or" sense (i.e. only one 
algorithm detecting a clear three jet structure is suffi- 
cient). 

The inbalance expressed by the three momentum px 
sum scales with the px scale the jet algorithm recon- 
structs to. Therefore the comparison is made against the 
k± algorithm with the momentum sum convoluted by 
the PTfiitei/pT,k ± distribution function, i.e. if one would 
"simulate" the filter behavior by apply the distribution 
function event by event to the kj_ algorithm. 

While not explicitly shown, we also checked our results 
against herwig 5.6.10 [18j |. and found no algorithm spe- 
cific, dissimilar behaviors, beside differences among the 
event generators that are reproduced by all algorithms. 

IV. DISCUSSION 

The prominent feature of jet reconstruction by filter- 
ing is the effect of the angular weight. It has the effect 
that the reconstructed jet px is shifted for wide angle 
fragmentation, while becoming less noticeable for sharply 
focused, high pt jets (Figured]). This accounts for the 
performance of the filter algorithm in the rejection of the 
soft background, which usually only contributes signifi- 
cantly to the event px when summed over large angles. 
This inherent discrimination against the background pro- 
vide us with a good start position to apply the formalism 
presented here to stronger background levels, including 
those usually found in heavy ion collisions. And while 
this effect complicates the determination of jet px, the 
parameter dependent shift of jet px and the subsequent 
need for calibration is present with all jet reconstruction 
algorithms. 

To show the insensitivity to the background, we cal- 
culated both the jet multiplicity and jet px specturm at 
8GeV/c initial hard scattering. In term of the jet mul- 
tiplicity (Figure [3]) , the filtering reproduces more accu- 
rately an initial hard scattering like event shape, which 
is dominated by two and three jet events. 

The px spectra (Figured]) shows a two component mix- 
ture, consisting of a high jet cross section, low px expo- 
nential distribution resulting from the soft QCD back- 
ground, and a peak from the triggered, initial hard scat- 
tering. Here both the k± algorithm and SISCone in 
fact produce a large amount of "jets" from background 



particles, thus creating the large jet multiplicity visi- 
ble in Figure [3] while the Gaussian filter suppresses the 
background approximately by a factor of 2 at an initial 
hard scattering of 8 GeV/c, thus producing a significantly 
more isolated and prominent peak for the trigger. This 
effect is consistent over a range of filter sizes and PT,trig? 
and becomes even stronger with increasing filter size or 
at lower p-r.trig- While SISCone e.g. fails to resolve the 
hard scattering cross section peak at 4GeV/c from the 
background for any choice of R < \ , this is possible with 
the filter algorithm and a > 0.8, albeit not at a = 0.5. 

In di-jet events, the Gaussian filter reproduces the 
angular behavior that previous jet reconstruction algo- 
rithms exhibit (Figures [5]) , and a slight difference in the 
three momentum sum is visible, with the filter having a 
slightly weaker exponential tail (Figure . The similari- 
ties are expected, since with adequate jet reconstruction 
performance, these variables should be dominated by the 
physical processes rather than the jet definitions. There 
is also a slight difference in the three jet opening angle 
when compared to the k± algorithm (J7J, and difference 
to the k± algorithm is smaller than the k± algorithm to 
SISCone. 



V. CONCLUSION 

In this paper, we presented a new algorithm for jet 
reconstruction by filtering, that beside essential pQCD 
properties of collinear and infrared safety, provides ro- 
bustness against soft background. Applied on event gen- 
erator final state particles, we showed that its perfor- 
mance in key areas matches the k± algorithm and SIS- 
Cone, such as for the di-jet balance and three jet open- 
ing angle, which are given by the physical jet production 
and jet fragmentation and not significantly modified by 
the jet definition. While the px scale differs from the 
reconstructed scale of the k± and cone algorithms, we 
note that an energy calibration is usually required for 
all algorithms for the extraction of the accurate parton 
energy. We also demonstrated, in this paper for p+p col- 
lisions, its superior ability in discriminating jets from the 
background, thus rejecting the associated, spurious jets 
contaminating events at low px- We hope to apply the 
algorithm presented here to most of the collision systems 
available at RfflC and LHC, thus providing an unified 
jet definition suitable across a large range of species, col- 
lision energies and jet px that so far was not attainable 
with the traditional jet reconstruction algorithms. 
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